////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// This file is a part of PMCG project. PMCG is the acronym of Physical       //
// Methods in Computer Graphics.                                              //
//                                                                            //
// Classical algorithms in computer graphics use discrete objects in 3D space //
// like vectors, matrices, triangles, planes. Physical way use an  analytical //
// objects like function, scalar and vector fields. Physical way gives new    //
// look to classical algorithms and allows to overcome some difficulties of   //
// pure geometric methods.                                                    //
//                                                                            //
// Project site:                                                              //
//     http://code.google.com/p/pmcg/                                         //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

#include <cmath>
#include "GeometryAlgs.h"

bool IsPointInPolyhedron(const TriMesh & mesh, const Vector3D & point)
{
    const double cFluxEps = 0.1;
    const double cDistanceEps = 0.002;

    double surfFlux = 0.0;

    TriMesh::Faces::const_iterator it;
    for (it = mesh.m_faces.begin(); it != mesh.m_faces.end(); ++it)
    {
        // face vertices
        const Vector3D & v0 = mesh.m_vertices[(*it).m_i0];
        const Vector3D & v1 = mesh.m_vertices[(*it).m_i1];
        const Vector3D & v2 = mesh.m_vertices[(*it).m_i2];

        // center of the face
        Vector3D c = (v0 + v1 + v2) / 3;

        // area of the face and outward normal
        Vector3D v10 = v1 - v0;
        Vector3D v20 = v2 - v0;
        Vector3D n = v10 ^ v20;
        double ds = sqrt(n * n) / 2;
        n.Normalize();

        // calc flux through the face
        Vector3D r = c - point;
        double l = sqrt(r * r);
        if (l < cDistanceEps)
            return true;
        double flux = (r * n) / (l * l * l) * ds;

        surfFlux += flux;
    }
    surfFlux /= (4 * M_PI);

    if (fabs(surfFlux - 1) < cFluxEps)
        return true;
    else
        return false;
}
